#include "common.h"

/*
 * GPU kernel to calculate sum and sum of square 
 */
__global__ void sum_sum2_kernel(float *out_sum, float *out_sum2, float *in, int n, int m) {
   __shared__ float sSumPtr[TILESIZE][TILESIZE];
   __shared__ float sSum2Ptr[TILESIZE][TILESIZE];
   int i;
   int offset;
   float x;

   i = blockIdx.y*blockDim.y + threadIdx.y; 
   if (i >= n) {
      return;
   }

   sSumPtr[threadIdx.y][threadIdx.x] = 0;
   sSum2Ptr[threadIdx.y][threadIdx.x] = 0;
  
   for (offset = 0; offset < m; offset += blockDim.x) {
      if (threadIdx.x + offset < m) {
         x = in[(blockIdx.y*blockDim.y + threadIdx.y)*m + threadIdx.x + offset];
         sSumPtr[threadIdx.y][threadIdx.x] += x;
         sSum2Ptr[threadIdx.y][threadIdx.x] += x*x;
      }
   }
   __syncthreads();

   for (unsigned int stride = TILESIZE/2; stride > 0; stride = stride >> 1) {
      if (threadIdx.x < stride) {
         sSumPtr[threadIdx.y][threadIdx.x] += sSumPtr[threadIdx.y][threadIdx.x+stride];   
         sSum2Ptr[threadIdx.y][threadIdx.x] += sSum2Ptr[threadIdx.y][threadIdx.x+stride];   
      }
      __syncthreads();
   }

   if (threadIdx.x == 0) {
      out_sum[i] = sSumPtr[threadIdx.y][0];
      out_sum2[i] = sSum2Ptr[threadIdx.y][0];
   }
}

/*
 * GPU kernel to subtract vector B from A 
 * Work only on 1D grid
 */
__global__ void sub_kernel(float *A, float *B, int size) {
   int j;

   j = blockIdx.x*blockDim.x + threadIdx.x;

   if (j < size) {
      A[j] -= B[j];
   }
}



